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Abstract 

New efficient and accurate numerical methods are proposed to compute ground 
states and dynamics of dipolar Bose-Einstein condensates (BECs) described by a 
three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a dipolar interaction 
potential. Due to the high singularity in the dipolar interaction potential, it brings 
significant difficulties in mathematical analysis and numerical simulations of dipolar 
BECs. In this paper, by decoupling the two-body dipolar interaction potential into 
short-range (or local) and long-range interactions (or repulsive and attractive interac- 
tions) , the GPE for dipolar BECs is reformulated as a Gross-Pitaevskii-Poisson type 
system. Based on this new mathematical formulation, we prove rigorously existence 
and uniqueness as well as nonexistence of the ground states, and discuss the existence 
of global weak solution and finite time blowup of the dynamics in different parameter 
regimes of dipolar BECs. In addition, a backward Euler sine pseudospectral method 
is presented for computing the ground states and a time-splitting sine pseudospectral 
method is proposed for computing the dynamics of dipolar BECs. Due to the adaption 
of new mathematical formulation, our new numerical methods avoid evaluating inte- 
grals with high singularity and thus they are more efficient and accurate than those 
numerical methods currently used in the literatures for solving the problem. Extensive 
numerical examples in 3D are reported to demonstrate the efficiency and accuracy of 
our new numerical methods for computing the ground states and dynamics of dipolar 
BECs. 

Key Words: Dipolar Bose-Einstein condensate, Gross-Pitaevskii equation, dipolar in- 
teraction potential, Gross-Pitaevskii-Poisson type system, ground state, backward Euler 
sine pseudospectral method, time-splitting sine pseudospectral method. 

1 Introduction 

Since 1995, the Bose-Einstein condensation (BEG) of ultracold atomic and molecular gases 
has attracted considerable interests both theoretically and experimentally. These trapped 
quantum gases are very dilute and most of their properties are governed by the interactions 
between particles in the condensate [31]. In the last several years, there has been a quest 
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for realizing a novel kind of quantum gases with the dipolar interaction, acting between 
particles having a permanent magnetic or electric dipole moment. A major breakthrough 
has been very recently performed at Stuttgart University, where a BEC of ^^Cr atoms has 
been realized in experiment and it allows the experimental investigations of the unique 
properties of dipolar quantum gases [22]. In addition, recent experimental developments 
on cooling and trapping of molecules [17], on photoassociation [43], and on Feshbach 
resonances of binary mixtures open much more exciting perspectives towards a degenerate 
quantum gas of polar molecules [35]. These success of experiments have spurred great 
excitement in the atomic physics community and renewed interests in studying the ground 
states [36, 48, 20, 21, 23, 34] and dynamics [25, 30, 32, 50] of dipolar BECs. 

At temperature T much smaller than the critical temperature Tc, a dipolar BEC is 
well described by the macroscopic wave function = t) whose evolution is governed 
by the three-dimensional (3D) Gross-Pitaevskii equation (GPE) [48, 36] 



— + y(x) + c/oivi' + (T^dip * IV'I') 



ip, xG M^ t > 0, (1.1) 



where t is time, x = (x, y, z) G M is the Cartesian coordinates, % is the Planck constant, 
m is the mass of a dipolar particle and V^(x) is an external trapping potential. When a 
harmonic trap potential is considered, F(x) = '^{uJ^x^ + uiytj^ + uJlz^) with uj^., Uy and uoz 

being the trap frequencies in x-, y- and z-directions, respectively. Uq = ^^^"^ describes 
local (or short-range) interaction between dipolcs in the condensate with the s-wave 
scattering length (positive for repulsive interaction and negative for attractive interaction). 
The long-range dipolar interaction potential between two dipoles is given by 



V^dip(x) 



^o/^dip 1 - 3(x • n 



|2/|X|2 



^ofidip l-3cos2(^) 



X ' 



X G 



X ■ 



(1.2) 



where /xq is the vacuum magnetic permeability, //dip is permanent magnetic dipole moment 
(e.g. Hdip = for ^^Cr with /x^ being the Bohr magneton), n = (ni, n2-,n^)^ G is the 

dipole axis (or dipole moment) which is a given unit vector, i.e. jnj = ^Jn\ + n'2+n^ = 1, 
and 9 is the angle between the dipole axis n and the vector x. The wave function is 
normalized according to 

(1.3) 



with Wo = min{a;a;,a;j^,a;^}, x 
, we obtain the dimensionless GPE in 3D from (1.1) as 



|V(x,t)|2 dx = 7V, 

where N is the total number of dipolar particles in the dipolar BEC. 
By introducing the dimensionless variables, t 

aox with oo = 

[48, 49, 31, 5]: 

idftpix, t) = 



muio ' ^ 



3/2 



- V + F(x) + + A (c/dip * H') 



X G 



t > 0, (1.4) 



where /3 



NUo 



A 



4(7^a;^+7^y^-|-7f z^) is the dimensionless 



harmonic trapping potential with 7a; = 7j/ = ^ and 7^ = and the dimensionless 



wo 



long-range dipolar interaction potential i7dip(x) is given as 

3 l-3(x -11)2/1x12 3 l-3cos2(0) 



^^dip(x) 



47r 



X ' 



47r 



X G 



X ' 



(1.5) 
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From now on, we will treat (3 and A as two dimensionless real parameters. We understand 
that it may not physical meaningful when A < for modeling dipolar BEC. However, it is 
an interesting problem to consider the case when A < at least in mathematics and it may 
make sense for modeling other physical system. In fact, the above nondimensionlization 
is obtained by adopting a unit system where the units for length, time and energy are 
given by ao, I/wq and Tiooq^ respectively. Two important invariants of (1.4) are the mass 
(or normalization) of the wave function 



iV(V^(-,t)):=||V(-,i)f = / |V(x,i)|'(ix= / |V^(x,0)|2dx=l, t>0 
and the energy per particle 



(1.6) 



= E{,P{;0)), t>0. 



dx 



(1.7) 



To find the stationary states including ground and excited states of a dipolar BEC, we 
take the ansatz 



V'(x,i) 



X e 



t > 0, 



(1.8) 



where /i G M is the chemical potential and (j) := (/>(x) is a time-independent function. 
Plugging (1.8) into (1.4), we get the time-independent GPE (or a nonlinear eigenvalue 
problem) 



/i(?!)(x) 



_lv2 + F(x) + ^|<^|2 + A (c/dip * l-^l') 



0(x), 



X G 



under the constraint 



(x)l^ dx 



(1.9) 



:i.io) 



The ground state of a dipolar BEC is usually defined as the minimizer of the following 
nonconvex minimization problem: 
Find (t)g e S and fj,^ eR such that 



min E( 



a ■- 



E3 := E{4>g 

where the nonconvex set S is defined as 

S:={ct>\ Uf = l,E{ct>)<<^] 
and the chemical potential (or eigenvalue of (1.9)) is defined as 



(1.11) 



:i.i2) 



|V0|2 + y(x)|<^|2 + p\ct>\^ + A (c/dip * 



dx 



= £;(^) + ijrj^|<^|4 + A(c/dip*|0l') 



dx. 



(1.13) 



In fact, the nonlinear eigenvalue problem (1.9) under the constraint (1-10) can be viewed 
as the Euler-Lagrangian equation of the nonconvex minimization problem (1.11). Any 
eigenfunction of the nonlinear eigenvalue problem (1.9) under the constraint (1.10) whose 
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energy is larger than that of the ground state is usuaUy called as an excited state in the 
physics literatures. 

The theoretical study of dipolar BECs including ground states and dynamics as well 
as quantized vortices has been carried out in recent years based on the GPE (1.1). For the 
study in physics, we refer to [16, 18, 24, 33, 1, 19, 24, 27, 28, 44, 45, 49, 52] and references 
therein. For the study in mathematics, existence and uniqueness as well as the possible 
blow-up of solutions were studied in [12], and existence of solitary waves was proven 
in [2]. In most of the numerical methods used in the literatures for theoretically and/or 
numerically studying the ground states and dynamics of dipolar BECs, the way to deal with 
the convolution in (1.4) is usually to use the Fourier transform [25, 20, 34, 46, 10, 41, 51]. 
However, due to the high singularity in the dipolar interaction potential (1.5), there are 
two drawbacks in these numerical methods: (i) the Fourier transforms of the dipolar 
interaction potential (1.5) and the density function are usually carried out in the 
continuous level on the whole space (see (2.5) for details) and in the discrete level on 
a bounded computational domain respectively, and due to this mismatch, there is a 
locking phenomena in practical computation as observed in [34] ; (ii) the second term in the 
Fourier transform of the dipolar interaction potential is g-type for 0-mode, i.e when ^ = 
(see (2.5) for details), and it is artificially omitted when ^ = in practical computation 
[34, 21, 29, 50, 49, 46, 10] thus this may cause some numerical problems too. The main 
aim of this paper is to propose new numerical methods for computing ground states and 
dynamics of dipolar BECs which can avoid the above two drawbacks and thus they are 
more accurate than those currently used in the literatures. The key step is to decouple the 
dipolar interaction potential into a short-range and a long-range interaction (sec (2.4) for 
details) and thus we can reformulate the GPE (1.4) into a Gross-Pitaevskii-Poisson type 
system. In addition, based on the new mathematical formulation, we can prove existence 
and uniqueness as well as nonexistence of the ground states and discuss mathematically 
the dynamical properties of dipolar BECs in different parameter regimes. 

The paper is organized as follows. In section 2, we reformulate the GPE for a dipolar 
BEG into a Gross-Pitaevskii-Poisson type system and study analytically the ground states 
and dynamics of dipolar BECs. In section 3, a backward Euler sine pseudospectral method 
is proposed for computing ground states of dipolar BECs; and in section 4, a time-splitting 
sine pseudospectral (TSSP) method is presented for computing the dynamics. Extensive 
numerical results are reported in section 5 to demonstrate the efficiency and accuracy of 
our new numerical methods. Finally, some conclusions are drawn in section 6. Throughout 
this paper, we adapt the standard Sobolev spaces and their corresponding norms. 

2 Analytical results for ground sates and dynamics 

Let r = |x| = ^Jx"^ + y'^ + and denote 



5n = n • V = nidx + n2dy + nsd- 



5, 



'nn — 



5n(9n). 



(2.1) 



Using the equality (see [30] and a mathematical proof in the Appendix) 





(2.2) 
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with (5(x) being the Dirac distribution function and introducing a new function 



we obtain 



C^dip * \^{;t)f = -|V'(x,i)|^ - 35„„ ((^(x,t)) , X G 



i > 0. 



(2.4) 



In fact, the above equaUty decouples the dipolar interaction potential into a short-range 

and a long-range interaction which correspond to the first and second terms in the right 
hand side of (2.4), respectively. In fact, from (2.1)-(2.4), it is straightforward to get the 
Fourier transform of [/dip(x) as 



(i7dip)(0 = -i + 



3(n-e)^ 



(2.5) 



Plugging (2.4) into (1.4) and noticing (2.3), we can reformulate the GPE (1.4) into a 
Gross-Pitaevskii-Poisson type system 



^V' + F(x) + XM{yi,t)f - 3Aa„„(^(x,i) 



VV(x,i) = -|V'(x,i)P, limv'(x,t)=0 x G M^, 

|x|— >-oo 



V'(x,i), (2.6) 
t > 0. (2.7) 



Note that the far-field condition in (2.7) makes the Poisson equation uniquely solvable. 
Using (2.7) and integration by parts, we can reformulate the energy functional E(-) in 
(1.7) as 



1 1 SA 



dx , 



(2.8) 



where tp is defined through (2.7). This immediately shows that the decoupled short- 
range and long-range interactions of the dipolar interaction potential are attractive and 
repulsive, respectively, when A > 0; and arc repulsive and attractive, respectively, when 
A < 0. Similarly, the nonlinear eigenvalue problem (1.9) can be reformulated as 



+ y(x) + (/3 - A) |</)|2 - 3Aa„n99(x) 



x 



VV(x) = -|<^(x)|2, xG 



lim <^(x) = 0. 

|x|-^oo 



(2.9) 
(2.10) 



2.1 Existence and uniqueness for ground states 

Under the new formulation for the energy functional E{-) in (2.8), we have 

Lemma 2.1 For the energy E{-) in (2.8), we have 

(i) For any (p E S, denote p(x) = |0(x)p for x G M^, then we have 

E{<l))>Em=E{^p), ycf>eS, 

so the minimizer (pg of (1.11) is of the form e*^"|0g| for some constant Oq G R. 
(a) When P > and —\P < A < /3, the energy E{y/p) is strictly convex in p. 



(2.11) 
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Proof: For any (j) ^ S, denote p = |0p and consider the Poisson equation 



VV(x) 



-I0(x) 



-P(x) 



X G 



lim (^(x) 
|x|-^oo 



0. 



Noticing (2.1) with |n| = 1, we have the estimate 

WdnV^h < ||1?V||2 = ||VV||2 = ||p||2 = 

(i) Write (/)(x) = e*^W|(/>(x)|, noticing (2.8) with V 



/ 



with D"^ = VV. 
) and (2.12), we get 
3A, 



V|<^| p + |<^nve(x)p + F(x)|</>|2 + -(/3 - A)|</>|^ + ^l^nV^r 



> 



VI 



+ y(x 



^ + ^(/3-A)M^ + Y|5nVv.p 
/0 G 5, 



(2.12) 
(2.13) 
c/x 

(2.14) 



and the equahty holds iff V0(x) = for x G M^, which means 0(x) = is a constant. 

(ii) From (2.8) with = cp and noticing (2.12), we can spht the energy E (y^) into 
two parts, i.e. 



E{^) = Eii^) + E2i^), 



where 



EiiVp)= I \\VVp? + V{x)p 



dx, 



^(/3-A)|p|^ + y|a„Vv.|^' 



dx. 



(2.15) 

(2.16) 

(2.17) 



As shown in [26], Ei {y/p) is convex (strictly) in p. Thus we need only prove E2 (y/p) is 
convex too. In order to do so, consider ^ypi G S, y/p2 £ S, and let ipi and (p2 be the 
solutions of the Poisson equation (2.12) with p = pi and p = p2, respectively. For any 
a G [0, 1], we have \Ja.p\ + (1 — a)p2 G -S, and 



aE2{y/p{) + (1 - a)E2{y/p^) - E2 \l api + (1 - a)p2 



= a{l — a) I 



1 3A 

-(^ - A)(pi - P2)' + Y|9nV(<^l - <^2)P 



dx, 



(2.18) 



which immediately implies that E2{y/p) is convex if /3 > and 0<A</5. If/5>0 and 
— < A < 0, noticing that aifi + (1 — a)ip2 is the solution of the Poisson equation (2.12) 
with p = api + {\ — a)p2, combining (2.13) with = (^1 — (^2 and (2.18), we obtain E2{-s/p) 
is convex again. Combining all the results above together, the conclusion follows. □ 

Now, we are able to prove the existence and uniqueness as well as nonexistence results 
for the ground state of a dipolar BEC in different parameter regimes. 

Theorem 2.1 Assume V{x) > for x G and lim V{x) = 00 (i.e., confining poten- 

|x|->oo 

tial), then we have: 

(i) If (3 > and < A < there exists a ground state (j)g G S, and the positive 
ground state \(f)g\ is unique. Moreover, (j)g = e'^^°\(f)g\ for some constant 9q G R. 

(ii) If P < Q, or /3 > and X < —^P or X > /3, there exists no ground state, i.e., 
inf £;(0) = -00. 
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Proof: (i) Assume /3 > and —^P < A < ^, we first show E{(j)) is nonnegative in S, i.e. 



1 Q \ 



dx>0, V(^g5. (2.19) 



In fact, when /3 > and < A < ^, noticing (2.8) with ip = (j), it is obvious that (2.19) is 
valid. When /3 > and < A < 0, combining (2.8) with ^ = (p, (2.12) and (2.13), we 
obtain (2.19) again as 



m > I 



|V(/.|2 + F(x 



dx 



|V<^P + y(x)|0|2 + -(;9 + 2A)|</.|4 



dx > 0. 



(2.20) 



Now, let {0"}^Q C S" be a minimizing sequence of the minimization problem (1.11). Then 
there exists a constant C such that 



iiv<^"ii2<c, \\r\u<c, 



Jm 



^dx <C, n > 0. 



(2.21) 



Therefore 0" belongs to a weakly compact set in L^, = {cp \ ||^||2 + ||V^||2 < oo}, 
and Ly = {(j) \ /j^a y(x)|(^(x)p dx < oo} with a weighted L^-norm given by ||i?!>||y = 
[/]R3 |</'(x)pT^(x)dx]^/^. Thus, there exists a (f)°° € //^ f] -^y Pi ^-^id a subsequence (which 
we denote as the original sequence for simplicity), such that 



(/)°°, inL^nL^nL^, 



V(^" ^ V<^°°, in 



(2.22) 



Also, we can suppose that is nonnegative, since we can replace them with 1^"! , which 
also minimize the functional E. Similar as in [26], we can obtain ||(?!'°°||2 = 1 due to the 
confining property of the potential y(x). So, (^°° G S*. Moreover, the L^-norm convergence 
of ^" and weak convergence in (2.22) would imply the strong convergence 0" — )■ G L^. 
Thus, employing Holder inequality and Sobolev inequality, we obtain 



n\2 



iOO\2| 



n||l/2 
l6 + 



n 



oo, 



(2.23) 



2<CM' 

< C2(||v</>"||f + \\v4>'^\\l^')\\r - (l^'^h ^ 0, 

which shows p" = {(f)")^ p°° = {(jf^)^ G L^. Since E^i^/p) in (2.17) is convex and 
lower semi-continuous in p, thus E2(4>°°) < lim E2{(j)'"). For £'i in (2.16), Ei((/)°°) < 

n— >oo 

lim £"1(0") because of the lower semi-continuity of the H^- and L?/-norm. Combining 
the results together, we know E{(/)°°) < lim E{(p"'), which proves that (/)°° is indeed a 
minimizer of the minimization problem (1.11). The uniqueness follows from the strictly 
convexity of E{yfp) as shown in Lemma 2.1. 

(ii) Assume ^ < 0, or /3 > and \ < —\j5 on \ > j5. Without loss of generality, we 
assume n = (0, 0, 1)-^ and choose the function 

2 ■ 



^£1 



,£2(X) 



1 



1 



(27r£i)V2 (27r£2)V4 



exp 



2ei 



exp 



z 
2i^. 



X G 



(2.24) 



with £1 and £2 two small positive parameters (in fact, for general n G satisfies |n| = 1, 
we can always choose 7^ ni G and 7^ n2 G such that {ni, n2, n} forms an 
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orthonormal basis of and do the change of variables x = [x, y, z)"^ to y = (x • ni, x • 
n2, X- n)-^ on the right hand side of (2.8), the following computation is still valid). Taking 
the standard Fourier transform at both sides of the Poisson equation 



-VV.i,.2(x) 
we get 



lim 99,,,,, (x) =0, (2.25) 
|x|— >-oo 



Iepcr.(e) = /Q;(0, ^e^'- 

Using the Plancherel formula and changing of variables, we obtain 

1 ......o 1 r /_ _,,x2 



(2.26) 



|(9nV(^, 



£1,£2 112 



(27r)2 



(27r)3 hs ICP 



(27r)3£iV£i (|^,|2 + |^2P)-ff + 16 
By the dominated convergence theorem, we get 
0, 



{KiiOfd^. £i,£2>0. (2.27) 



£2/£l ->■ +00, 



£1,£2 II2 



(2.28) 



When fixed ei-y/£2, the last integral in (2.27) is continuous in e2/£i > 0. Thus, for any 
a G (0,1), by adjusting e2/ei := > 0, we could have ||9nV9?£i,£2 Hi = c>;||^£i,£2 Ill- 
Substituting (2.24) into (2.16) and (2.17) with ^fp = 4>e'L,e2 under fixed £2/£i > 0, we get 



£1,£2 

Jr3 l 



2Wei,£2J 



- (/3-X + 3aX))\(pei,e2\ = — 

2 Jm.3 2 £l^/£2 



(2.29) 

(2.30) 



with some constants Ci, C2, C3 > independent of £1 and £2. Thus, if /S < 0, choose 
a = 1/3; if ;8 > and A < choose 1/3 - < a < 1; and if /? > and A > /3, 

choose < a < 4 (1 — as £1, £2 ^ 0+, we can get inf £^(0) = lim Ei{(f)si,e2) + 

^ ' <?!>eS £l,£2-^0+ ' 

£^2 (</>£!, £2) = —00, which implies that there exists no ground state of the minimization 
problem (1.11). □ 

By splitting the total energy E{-) in (2.8) into kinetic, potential, interaction and dipolar 
energies, i.e. 

E{cl>) = E^i^icj)) + E^ot{(t>) + E;^t{<t>) + £;dip(<^), (2.31) 



where 



\ I |V0(x)|2dx, E^M = I V(x)|0(x)|2dx, E,^,{4>) = ^ / |</.(x)|4dx, 

/ (c/dip * |</.h |</)(x)|2dx = ^ / |</,(x)p [3 (5„„v.)^ - |0(x)pl 
JR3 ^ ^ 2 Jk3 L J 

/ f-|(/)(x)|4 + 3|a„V¥^|'ldx, (2.32) 

7m3 L J 

with defined in (2.10), we have the following Viral identity: 
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Proposition 2.2 Suppose (pe is a stationary state of a dipolar BEC, i.e. an eigenfunction 
of the nonlinear eigenvalue problem (1.9) under the constraint (1.10), then we have 

2£^kin(0e) - 2£;trap('/'e) + 3^int(</'e) + ^E^iMe) = 0- (2-33) 

Proof: Follow the analogous proof for a BEC without dipolar interaction [31] and we 
omit the details here for brevity. □ 

2.2 Analytical results for dynamics 

The well-posedness of the Cauchy problem of (1.1) was discussed in [12] by analyzing the 
convolution kernel C/dip(x) with detailed Fourier transform. Under the new formulation 
(2.6)-(2.7), here we present a simpler proof for the well-posedness and show finite time 
blow-up for the Cauchy problem of a dipolar BEC in different parameter regimes. Denote 

X = [ue H^{m.^) I \\ufx = hWh + ||V«||i2 + / y(x)|?x(x)|2 dx < GO 

Theorem 2.3 (Well-posedness) Suppose the real-valued trap potential V{yi) G C°°(M^) 
such that F(x) > /or x G D°1/(x) G L°°(M3) for all a G N|] with \a\ > 2. 

For any initial data il){^,t = 0) = V'o(x) G X , there exists Tmax G (0, +00] such that the 
problem (2.6)- (2.7) has a unique maximal solution G C {[0,Tfnax), X) . It is maximal in 
the sense that i/Tmax < 00, then \\ip{-,t)\\x — >■ oo when t — )• T^^x. Moreover, the mass 
N{^p{-,t)) and energy E{ilj{-,t)) defined in (1.6) and (1-7), respectively, are conserved for 
t G [OjTmax)- Specifically, if P > and —^P < A < ^, the solution to (2.6)-(2.7) is global 
in time, i.e., Tmax = 00. 

Proof: For any (p & X, let (p be the solution of the Poisson equation (2.12), denote 
p = and define 

G{cP,^):=G{p) = \ j |0(x)|25„„(^(x)dx, 5(</,) = ^^^ = 0a„n9^, (2.34) 

where / denotes the conjugate of /. Noticing (2.13), it is easy to show that G(0) G 
CHX,E), 5(0) G C{X,LP) for some p G (6/5,2], and 

\\g{u)-g{v)\\Lv <C{\\u\\x + \\v\\x)\\u-v\\Lr, for some r G [2, 6), "iu^v e X. (2.35) 

Applying the standard Theorems 9.2.1, 4.12.1 and 5.7.1 in [13, 40] for the well-posedness 
of the nonlinear Schrodinger equation, we can obtain the results immediately. □ 

Theorem 2.4 (Finite time blow-up) If (3 < 0, or P > and X < — or A > /?, 
and assume V(k) satisfies 3V{x.) + x • W{x.) > for x G M^. For any initial data 
ip^x^t = 0) = ^o(x) E X to the problem (2.6)-(2.7), there exists finite time blow-up, i.e., 
Tmax < 00, if one of the following holds: 

(i) S(^o) < 0; 

(a) E{'ijjo) = and Im {J^s V'o(x) (x • VV'o(x)) dx) < 0; 
(Hi) E{i/jo) > and Im (j^s V'o(x) (x • VV'o(x)) dx) < -^/3E(^po}\\x'tpo\\L2 ; 
where iToa(f) denotes the imaginary part of f. 
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Proof : Define the variance 



av{t):=av{'^{-,t))= ( \^\^\il^{^,t)\^ = 5:,{t) + 5y{t) + 5,{t), t>0, (2.36) 



where 



aa{t) ■- cra{^j{-,t)) = a^|V'(x,t)p(ix, a = x,y,z. (2.37) 

For a = X, or y or z, differentiating (2.37) with respect to t, noticing (2.6) and (2.7), 
integrating by parts, we get 

7 /I 

-(T„(t) = -i [aV^(x, t)dail^{yi, t) - aV(x, t)daijj{yi, t)] dx, t > 0. (2.38) 

at J-R3 

Similarly, we have 

^aa{t) = I [2|a«V|' + + eAIV-pa^aSnn.^ - 2a|V'p5ay(x)l dx. (2.39) 

Noticing (2.7) and 

- / VV (X • VannV') d^= \ j \dr,V^\^ dx, 
yM3 2 Jm3 

summing (2.39) for a = y and z, using (2.36) and (1.7), we get 

^C7y(t) = 2 (iVV'l^ + ^(^ - A)|V^|^ + ^AlSnVV'l^ - |Vp(x • Vy(x))) dx 

= QE{il))- I |V^(x,t)|2-2 / |V'(x,t)|2(3y(x) + x- Vy(x)) dx 

< 6£;(V') = 6£;(V'o), t > 0. (2.40) 

Thus, 

(Tv{f) < ^E{Tjjo)t^ + a'y{0)t + CTy(O), t > 0, 

and the conclusion follows in the same manner as those in [40, 13] for the standard non- 
linear Schrodinger equation. □ 



3 A numerical method for computing ground states 

Based on the new mathematical formulation for the energy in (2.8), we will present an effi- 
cient and accurate backward Euler sine pseudospectral method for computing the ground 
states of a dipolar BEC. 

In practice, the whole space problem is usually truncated into a bounded computational 
domain Q = [a, b] x [c, d] x [e, /] with homogeneous Dirichlet boundary condition. Various 
numerical methods have been proposed in the literatures for computing the ground states 
of BEC (see [37, 15, 4, 3, 7, 14, 11] and references therein). One of the popular and efficient 
techniques for dealing with the constraint (1.10) is through the following construction 
[4, 8, 3]: Choose a time step At > and set tn = n At for n = 0, 1, . . . Applying the 
steepest decent method to the energy functional E{(j)) in (2.8) without the constraint 
(1.10), and then projecting the solution back to the unit sphere S at the end of each time 
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interval [tn,tn+i] in order to satisfy the constraint (1.10). This procedure leads to the 
function (^(x, t) is the solution of the following gradient flow with discrete normalization: 



-V2 - F(x) - (/3 - A)|(^(x,t)|2 + 3A5„„(^(x,t) 



(x,t), 



VV(x,t) = -|^(x,i)p, 



(/>(x, t 



n+1) 



X,t 



X G tfi <t < fn+l) 

X G Jl, n > 0, 



</.(x,t)| 

^(x,o) = ?;>o(x), 



(/9(x, i)| 



0, 



with 1 1 



t > 0, 
1; 



(3.1) 
(3.2) 

(3.3) 

(3.4) 
(3.5) 



where ^(x, t^) = lim ^(x, t). 



Let M, K and L be even positive integers and define the index sets 

TMKL = {{j,k,l) |i = l,2,...,M-l, A; = l,2,...,i^-1, / = 1, 2, . . . , L - 1}, 
7MKL = {{j,k,l) \j = 0,l,...,M, k = 0,l,...,K, l = 0,l,...,L}. 

Choose the spatial mesh sizes as hx = ^j^, hy = and hz = ^-j^ and define 



Xj := a + j hx, yk = c + k hy, zi = e + I hz 
Denote the space 



{j,k,l)eTZ 



MKL- 



Ymkl = span{$jfe;(x), (j, k, I) € Tmkl}, 



with 



^jfe^x) = sin (/xj(a; - afj sin (/x^(y - c)) sin (/xf (z - e)) , xeCl, [j, k, I) G Tmkl, 



h — 



4 



d — c 



{j, k, I) G Tmkl; 



f-e 

and Pmkl - Y = {(p e C{Q,) | (^(x)|xean = 0} — >■ Ymkl be the standard project operator 
[38], i.e. 

M-lK-lL-l 

{PmKLv){^) = X] X] X] ^pqs(x), X G G 

p=l q=l s=l 

with 

= v{:x) ^pqs(x) dx, (p, q, s) G Tmxl- 

Then a backward Eulcr sine spectral discretization for (3.1)-(3.5) reads: 
Find 0"+i(x) G iMi^L (i-e. ^^^(x) G Ymkl) and <y9''(x) G Ymkl such that 



(3.6) 



<^+(x)^^r(x) ^ 1 v2^+(^) _ p^^^ I ^ _ A)|</,"(x)|2 + 3Aa„„^"(x)] .^+(x)} ,(3.7) 



VV"(x) = -Fmxl ([-^"(x)!') , </>"+^(x) - 



'^+(x) 
ll<^+(x)||2' 



x € f2, n > 0; 



(3.8) 



where (/)°(x) = Pmj^l (</'o(x)) is given. 

The above discretization can be solved in phase space and it is not suitable in prac- 
tice due to the difficulty of computing the integrals in (3.6). We now present an effi- 
cient implementation by choosing ^°(x) as the interpolation of 0o(x) on the grid points 
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{{xj,yk,zi), {j,k,l) G Tmkl}^ i-e 4'°ixj,yk, zi) = Mxj^yk,zi) for {j,k,l) G T^kl^ and 
approximating the integrals in (3.6) by a quadrature rule on the grid points. Let 0^^,^ 
and ip'ji^i be the approximations of (l){xj,yk, zi,tn) and ip{xj,yk, zi,tn), respectively, which 
are the solution of (3.1)-(3.5); denote pji^^ = |<^"fcj^ and choose (p^f^i = (j)Q{xj,yk, zi) for 

(j, k, I) G T^j^i^. For n = 0, 1, . . ., a backward Euler sine pseduospectral discretization for 
(3.1)-(3.5) reads: 



- {^l^%ki = = Plki. = ^li^, ij,k,l) e Tmkl, (3.10) 

-^ofeV = ^mW = = - = = 0, (j, fc, /) e TStKL, (3.11) 
V'ow = VMki = ^loi = ^Iki = ^]ko = ^]kL = 0, (j, k, I) e T^KL, (3.12) 

where and d^^ are sine pseudospectral approximations of and dnn, respectively, 
defined as 

M-l K-1 L-1 / ■ \ / , \ / , 

ISTT 



L 



p=l q=l s=l ^ ^ \ / \ 

(5^n^")l,., = E E E (..p , ( , (9nn$p,.(x))|(,^.,^^,,,) , (i, k, I) G TmK. , (3 

p=i g=i s=i ^^P' vA*?/ yH'sj 
with {(t>'^)pqg {{p, q, s) G Tmkl) the discrete sine transform coefficients of the vector ^" as 

('^U = ]i^EEE^;--(f^)-(^^ (p,.,a)GrM.. (3.14) 

and the discrete /i-norm is defined as 



M-17V-1L-1 

u+\\l = h,hyh,Y: EE 14^/1'- 

j=i k=i 1=1 

Similar as those in [6], the linear system (3.9)-(3.12) can be iteratively solved in phase 
space very efficiently via discrete sine transform and we omitted the details here for brevity. 



4 A time-splitting sine pseudospectral method for dynamics 

Similarly, based on the new Gross-Pitaevskii-Poisson type system (2.6)-(2.7), we will 
present an efficient and accurate time-splitting sine pseudospectral (TSSP) method for 
computing the dynamics of a dipolar BEG. 

Again, in practice, the whole space problem is truncated into a bounded computational 
domain Q = [a, b] x [c, d] x [e, /] with homogeneous Dirichlet boundary condition. From 
time t = into time t = tn+i, the Gross-Pitaevskii-Poisson type system (2.6)-(2.7) is solved 
in two steps. One solves first 

^5^V(x,^) = -^VV(x,^), x G J^, V(x, t)^^^^ = 0, tn<t<tn+i, (4.15) 
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for the time step of length At, foUowed by solving 

z5tV(x,t) = [y(x) + {13- A)|V'(x,t)|2 - 3Aa„„y^(x,t)] V(x,t), (4.16) 
VVx,t) = -|V'(x,t)P, ^en, tn<t<tn+i; (4.17) 
<^(x,i)lxe9n = 0' V'(x,i)Le9n = 0' in < t < Wi; (4-18) 

for the same time step. Equation (4.15) will be discretized in space by sine pseudospectral 

method and integrated in time exactly [9]. For t G [^ni^n+i]) the equations (4.16)-(4.18) 
leave \ip\ and if invariant in t [5, 9] and therefore they collapses to 

iatV(x,t)= [y(x) + (/3-A)|V'(x,i„)|2-3Aa„„¥)(x,t„)] V(x,t), xefi, tn<t<tn+l, (4.19) 
VV(x,t„) = -|V(x,Ur, xea (4.20) 

Again, equation (4.20) will be discretized in space by sine pseudospectral method [9, 38] 
and the linear ODE (4.19) can be integrated in time exactly [5, 9]. 

Let V'jfe; and ip^j^i be the approximations of tp{xj,yk, zi,tn) and ip{xj,yk, zi,tn), re- 
spectively, which are the solution of (2.6)-(2.7); and choose ijjjj^i = ijjo{xj,yk, zi) for 

{j,k,l) e Tmkl- For n = 0, 1,..., a second-order TSSP method for solving (2.6)-(2.7) 
via the standard Strang splitting is [39, 5, 9] 

= E E (H,.- (^) sin (^) sin (i^ 

p=l q=l 5=1 \ / \ / \ 

^g) ^ ^-iA4n.....)+(.-A)i.-i-3A(a,„..>)|^.j ^a)^ ^.^^^^ ^ ^^^^ ^^21) 
= E E E (V^W- (^) - (^) - (x 

p=l q=l 5=1 \ / \ / \ 



where {ip"')pqr and (V'^^^)p5r iiP^I^^) ^ Tmkl) are the discrete sine transform coeffi- 
cients of the vectors -0" and respectively (defined similar as those in (3.14)); and 
(dnn^'"^^) .^^ can be computed as in (3.13) with p^^^ = pj'^j^ := iV'jHp for (j, k, I) G Tmkl- 

The above method is explicit, unconditionally stable, the memory cost is 0{MKL) 
and the computational cost per time step is O {MKL\il{MKL)). In fact, for the stability, 
we have 

Lemma 4.1 The TSSP method (4-21) is normalization conservation, i.e. 

M-lK-lL-l M-1K~1L~1 
j=l k=l 1=1 j=l k=l 1=1 

Proof: Follow the analogous proof in [5, 9] and we omit the details here for brevity. □ 

5 Numerical results 

In this section, we first compare our new methods and the standard method used in the 
literatures [49, 46, 41, 10] to evaluate numerically the dipolar energy and then report 
ground states and dynamics of dipolar BECs by using our new numerical methods. 
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5.1 Comparison for evaluating the dipolar energy 

Let 



X 



Then the dipolar energy £'dip(0) in (2.32) can be evaluated analytically as [42] 



dip I 



47rv27r 



1+2k^ 3«;^ arctan ( Vt"^ — 1 ) 

0, 

1+2k-^ 1.5k^ 1^ ( l+yr^ 

(1-^2)71^ Vi-VT^ 



K > 1, 
K = 1, 
K < 1, 



(5.1) 



(5.2) 



with K = y^- This provides a perfect example to test the efficiency of different numerical 
methods to deal with the dipolar potential. Based on our new formulation (2.32), the 
dipolar energy can be evaluated via discrete sine transform (DST) as 

j=i k=i 1=1 L 

where {d^n'P"')\jki is computed as in (3.13) with = \(l){xj,yk, zi)\^ for {j,k,l) G T^kl- 
In the literatures [49, 41, 46, 10], this dipolar energy is usually calculated via discrete 
Fourier transform (DFT) as 

Xh h h ^^-i-^-i^-i 
Edipi^) ~ ''J' ' E E E l^Pi^j^Vk^zi) 

j=i k=i 1=1 

where J- and J-^^ are the discrete Fourier and inverse Fourier transforms over the grid 
points {{xj,yk, zi), {j,k,l) G T^j^i}, respectively [46]. We take A = 247r, the bounded 
computational domain f2 = [—16, 16]^, M = K = L and thus h = hx = hy = hz = ^■ 
Table 1 lists the errors e := -Edip(</') — E^ip with E^^^ computed numerically via either 
(5.3) or (5.3) with mesh size h for three cases: 

• Case I. jx = 0.25 and 7^ = 1 which implies k = 2.0 and E^ipict)) = 0.0386708614; 

• Case II. 7a; = 72 = 1 which implies k = 1.0 and £^dip(?^) = 0; 

• Case III. 7j; = 2 and 7^ = 1 which implies k = \/a5 and Edip{(j)) = -0.1386449741. 



^ikl ((C/dip)(2/x^,2/x^,2/x^) • Jp,,(|</.p)) 





Case I 


Case II 


Case III 


DST DFT 


DST DFT 


DST DFT 


M = 32kh = 1 
M = 64kh = 0.5 
M = 128kh = 0.25 


2.756E-2 2.756E-2 

1.629E-3 1.614E-3 
1.243E-7 1.588E-5 


3.555E-18 1.279E-4 

9.154E-18 1.278E-4 
7.454E-17 1.278E-4 


0.1018 0.1020 

9.788E-5 2.269E-4 
6.406E-7 1.284E-4 



Table 1: Comparison for evaluating dipolar energy under different mesh sizes h. 

From Tab. 1 and our extensive numerical results not shown here for brevity, we can 
conclude that our new method via discrete sine transform based on a new formulation is 
much more accurate than that of the standard method via discrete Fourier transform in 
the literatures for evaluating the dipolar energy. 
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5.2 Ground states of dipolar BECs 



By using our new numerical method (3.9)-(3.12), here we report the ground states of 
a dipolar BEC (e.g., ^^Cr [30]) with different parameters and trapping potentials. In 
our computation and results, we always use the dimensionless quantities. We take M = 
K = L = 128, time step At = 0.01, dipolar direction n = (0,0,1)"^ and the bounded 
computational domain Q = [—8,8]^ for all cases except Q = [—16,16]^ for the cases 
= 1, 5, 10 and fl = [-20, 20]^ for the cases = 50, 100 in Table 2. The 

ground state (bn is reached numerically when 1 1 6'^'^^ — (t>'^\\oo '■= max I (p'li, ^ — 

< e := 10-6 in (3.9)-(3.12). Table 2 shows the energy ■= E{(t)g), chemical 
potential jjP := n{(t>g), kinetic energy E^^^ := ^^kinC^s), potential energy E^^^ := Epotifpg), 
interaction energy Ef^^^ := E\a\.{4>g)-, dipolar energy E'^jp := E^^p{(f)g), condensate widths 
(t| := (Tx{(t>g) and erf := crzifpg) iii (2.37) and central density Pg{0) := |^g(0, 0, 0)p with 
harmonic potential V{x, y,z) = ^ {x^ + + 0.25z^) for different /3 = 0.20716A/' and A = 
0.033146A?^ with N the total number of particles in the condensate; and Table 3 lists 
similar results with /? = 207.16 for different values of —0.5 < ^ < 1. In addition. Figure 
1 depicts the ground state 0g(x), e.g. surface plots of \(f)g{x,Q, z)\'^ and isosurface plots of 
|^g(x)| = 0.01, of a dipolar BEC with (3 = 401.432 and A = 0.16/3 for harmonic potential 
F(x) = i (x^ + y2 + z'^Y double-well potential F(x) = \{x'^ + + z^) + 46-^^/^ and 
optical lattice potential y(x) = i (x^ + + z^) + 100 [sin^ (f .t) + sin^ (f y) + sin^ (f z)] ; 
and Figure ?? depicts the ground state <^g(x), e.g. isosurface plots of |(^g(x)| = 0.08, of 
a dipolar BEC with the harmonic potential V^(x) = ^ (x^ + J/^ + z^) and (5 = 207.16 for 
different values of— 0.5 < ^ < 1. 



N 
10000 


E9 






E^ 

-^pot 


E9, 

mt 






^ z 


Pg{^) 


0.1 


1.567 


1.813 


0.477 


0.844 


0.262 


-0.015 


0.796 


1.299 


0.06139 


0.5 


2.225 


2.837 


0.349 


1.264 


0.659 


-0.047 


0.940 


1.745 


0.02675 


1 


2.728 


3.583 


0.296 


1.577 


0.925 


-0.070 


1.035 


2.009 


0.01779 


5 


4.745 


6.488 


0.195 


2.806 


1.894 


-0.151 


1.354 


2.790 


0.00673 


10 


6.147 


8.479 


0.161 


3.654 


2.536 


-0.204 


1.538 


3.212 


0.00442 


50 


11.47 


15.98 


0.101 


6.853 


4.909 


-0.398 


2.095 


4.441 


0.00168 


100 


15.07 


21.04 


0.082 


9.017 


6.498 


-0.526 


2.400 


5.103 


0.00111 



Table 2: Different quantities of the ground states of a dipolar BEC for (3 = 0.20716A'^ and 
A = 0.033146A/' with different number of particles A^. 

From Tabs. 2&:3 and Figs. 1&:??, we can draw the following conclusions: (i) For fixed 
trapping potential F(x) and dipolar direction n = (0,0, 1)^, when (3 and A increase with 
the ratio ^ fixed, the energy E^ , chemical potential /x^, potential energy E^^^^, interaction 
energy Ef^^^, condensate widths (t| and erf of the ground states increase; and rcsp., the 
kinetic energy E^^^^ dipolar energy -E'jjp and central density Pg(0) decrease (cf. Tab. 2). 
(ii) For fixed trapping potential V^(x), dipolar direction n = (0, 0, 1)"^ and /3, when the ratio 
^ increases from —0.5 to 1, the kinetic energy -E^jn, interaction energy E?^^, condensate 
widths erf and central density Pg{Q) of the ground states increase; and resp., the energy £'^, 
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A 

? 








^pot 


TP3 


^dip 




rr9 




-0.5 


2.957 


3.927 


0.265 


1.721 


0.839 


0.131 


1.153 


1.770 


0.01575 


-0.25 


2.883 


3.817 


0.274 


1.675 


0.853 


0.081 


1.111 


1.879 


0.01605 





2.794 


3.684 


0.286 


1.618 


0.890 


0.000 


1.066 


1.962 


0.01693 


0.25 


2.689 


3.525 


0.303 


1.550 


0.950 


-0.114 


1.017 


2.030 


0.01842 


0.5 


2.563 


3.332 


0.327 


1.468 


1.047 


-0.278 


0.960 


2.089 


0.02087 


0.75 


2.406 


3.084 


0.364 


1.363 


1.212 


-0.534 


0.889 


2.141 


0.02536 


1.0 


2.193 


2.726 


0.443 


1.217 


1.575 


-1.041 


0.786 


2.189 


0.03630 



Table 3: Different quantities of the ground states of a dipolar BEC with different values 
of ^ with ^ = 207.16. 

chemical potential /i^, potential energy -Bpot, dipolar energy £^^jp and condensate widths 
a% decrease (cf. Tab. 3). (iii) Our new numerical method can compute the ground states 
accurately and efficiently (cf. Figs. 1&:??). 

5.3 Dynamics of dipoleir BECs 

Similarly, by using our new numerical method (4.21), here we report the dynamics of 
a dipolar BEC (e.g., ^^Cr [30]) under different setups. Again, in our computation and 
results, we always use the dimensionless quantities. We take the bounded computational 
domain Q = [-8,8]^ x [-4,4], M = K = L = 128, i.e. h = h^c = hy = 1/8, = 1/16, 
time step At = 0.001. The initial data '0(x, 0) = V'o(x) is chosen as the ground state 
of a dipolar BEC computed numerically by our numerical method with n = (0,0,1)-^, 
V{x) = i(a;2 y2 + 25z^), (3 = 103.58 and A = 0.8/3 = 82.864. 

The first case to study numerically is the dynamics of suddenly changing the dipolar 
direction from n = (0,0,1)-^ to n = (1,0,0)-^ at t = and keeping all other quan- 
tities unchanged. Figure ?? depicts time evolution of the energy E{t) := E(^p(-,t)), 
chemical potential fi{t) = ^{ilj{-,t), kinetic energy -Ekin(i) '■= E]^m{ipi',t)), potential 
energy Epot{t) := Epot{i'{-,t)), interaction energy Eint{t) := EiT^t{'ip{-,t)), dipolar en- 
ergy Edip{t) := Edip{t/j{-,t)), condensate widths ax{t) := ax{il^{-,t)), az{t) := az{tp{-,t)), 
and central density p{t) := \ip{0,t)\'^, as well as the isosurface of the density function 
/9(x, t) := ['(/'(x, t)p = 0.01 for different times. In addition. Figure ?? show similar results 
for the case of suddenly changing the trapping potential from V^(x) = ^(a;^ + + 25^:^) to 
y(x) = ^{x'^+y'^+^z'^) at t = 0, i.e. decreasing the trapping frequency in z-direction from 
5 to |, and keeping all other quantities unchanged; Figure ?? show the results for the case 
of suddenly changing the dipolar interaction from A = 0.8/3 = 82.864 to A = 4/3 = 414.32 
at t = while keeping all other quantities unchanged, i.e. collapse of a dipolar BEC; and 
Figure ?? show the results for the case of suddenly changing the interaction constant /3 
from P = 103.58 to P = —569.69 at t = while keeping all other quantities unchanged, 
i.e. another collapse of a dipolar BEC. 

From Figs. ??, ??, ?? and ??, we can conclude that the dynamics of dipolar BEC can 
be very interesting and complicated. In fact, global existence of the solution is observed 
in the first two cases (cf. Figs. ??&??) and finite time blow-up is observed in the last 
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Figure 1: Surface plots of \cl)g{x,0, z)\'^ (left column) and isosurface plots of \(j)g{x,y, z)\ = 
0.01 (right column) for the ground state of a dipolar BEC with /3 = 401.432 and A = 0.16/3 
for harmonic potential (top row), double- well potential (middle row) and optical lattice 
potential (bottom row). 

two cases (cf. Figs. ??&??). The total energy is numerically conserved very well in our 
computation when there is no blow-up (cf. Figs. ??&??) and before blow-up happens (cf. 
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Figs. ??&??). Of course, it is not conserved numerically near or after blow-up happens 
because the mesh size and time step are fixed which cannot resolve the solution. In 
addition, our new numerical method can compute the dynamics of dipolar BEC accurately 
and efficiently. 

Remark 5.1 Due to size limit at ariv, to read the full figures, you can download this 
paper from: 

http://ww'w. math. nus. edu.sg/~bao/PS/dipolar-bec.pdf 

6 Conclusions 

Efficient and accurate numerical methods were proposed for computing ground states 
and dynamics of dipolar Bose-Einstein condensates based on the three-dimensional Gross- 
Pitaevskii equation (GPE) with a nonlocal dipolar interaction potential. By decoupling 
the dipolar interaction potential into a short-range and a long-range part, the GPE for a 
dipolar BEC is rc-formulatcd to a Gross-Pitaevskii-Poisson type system. Based on this 
new mathematical formulation, we proved rigorously the existence and uniqueness as well 
as nonexistence of the ground states, and discussed the dynamical properties of dipolar 
BEC in different parameter regimes. In addition, the backward Euler sine pseudospectral 
method and time-splitting sine pseudospectral method were proposed for computing the 
ground states and dynamics of a dipolar BEC, respectively. Our new numerical methods 
avoided taking the Fourier transform of the nonlocal dipolar interaction potential which is 
highly singular and causes some numerical difficulties in practical computation. Compari- 
son between our new numerical methods and existing numerical methods in the literatures 
showed that our numerical methods perform better. Applications of our new numerical 
methods for computing the ground states and dynamics of dipolar BECs were reported. 
In the future, we will use our new numerical methods to simulate the ground states and 
dynamics of dipolar BEC with experimental relevant setups and extend our methods for 
rotating dipolar BECs. 
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Appendix Proof of the equality (2.2) 

Let 
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For any n G satisfies |n| = 1, in order to prove (2.2) holds in the distribution sense, it 
is equivalent to prove the following: 



/ </,(x)/(x)dx = -^/(0)- / /(x) a„„ (-) dx, V/(x) e Co~(K3). (A.2) 
Jm3 o Jm.3 \rj 

I |x| > e}. It is 



For any fixed e > 0, let = {x G | |x| < e} and i?^ = {x G 
straightforward to check that 



^(x) 



Ot^xG 



(A.3) 



Using integration by parts and noticing (A.3), we get 
/ <^(x)/(x)dx = - / /(x) 5„„ (-) dx 

- k (^) '"^^^^^^ ^ L ^^^^ ^ (^) 

= - [ -5„n(/(x))o!x + /f + /|, 



(A.4) 



where 



n • X 



/(x) 



a„ /| := - ^ a„ (/(x)) (A.5) 



Prom (A.5), changing of variables, we get 

^2 



n • X 



as. 



-/(x)dS = 



(n • x)^ 



/(ex) e^dS 



dBi 



(n • x)V(O) d5 - / (n • x)^ [/(ex) - /(O)] 



9Bi 



(A.6) 



Choosing 7^ ni G and 7^ n2 G such that {ni, n2 n} forms an orthornormal basis 
of R^, by symmetry, we obtain 



JdBi o Jd 



3 JdBi 



3 JdBi 

dS 



(n • x)^ + (ni • x)2 + (na • x)^ 

47r 



3 JdBi 



I (n • x)2 (/(ex) - /(O)) d5 = / (n • ^fe [x • V/(^ex)] 

JaSi JdBi 

< e ||V/||ioo(s,) / dS< 47re ||V/|Uo.(s,), 



where < < 1. Plugging (A.7) and (A. 8) into (A.6), we have 



dS 



dS 



(A.7) 



(A. 



(A.9) 
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Similarly, for e — )• 0"*" , we get 




dx < 27r£2 ||L>V||loc(b,) ^ 0, 



(A.10) 



(A.ll) 



Combining (A.9), (A. 10) and (A.ll), taking e -)• 0+ in (A.4), we obtain 



/ </.(x)/(x)dx 



fm)-f 

o Jm 



R3 r 



- 5n„(/(x)) dx, V/(X) G C^{R^). 



(A.12) 



Thus (A. 2) follows from (A.12) and the definition of the derivative in the distribution 



and the equality (2.2) is proven. □ 
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